Setting

path_proj = here::here()
path_source = file.path(path_proj, "source")

source(file.path(path_source, "simulation", "simulations_functions_final.R"))
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:reshape2':
## 
##     smiths
source(file.path(path_source, "functions", "plot_function.R"))
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
source(file.path(path_source, "functions", "fit_function.R"))
source(file.path(path_source, "functions", "table_function.R"))

# place for draws
# mac
posterior_draws_path = file.path(Sys.getenv("HOME"), "Desktop", "draws", "testEach")
# windows
#posterior_draws_path = file.path(Sys.getenv("USERPROFILE"), "Desktop", "draws", "testEach")

#data path
data_save_path = file.path(path_proj, "data", "fitted_model", "simulation", "1. q_constant")
#models
q_constant <- file.path(path_proj, "source", "models", 
                      "q-constant.stan")
b_constant <- file.path(path_proj, "source", "models", 
                     "b-constant.stan")
b_rw <- file.path(path_proj, "source", "models", 
                     "b-rw1.stan")
b_ou <-  file.path(path_proj, "source", "models", 
                      "b-ou.stan")

compiled_models <- list(
  q_constant = cmdstan_model(q_constant),
  b_constant = cmdstan_model(b_constant),
  b_rw = cmdstan_model(b_rw),
  b_ou = cmdstan_model(b_ou)
)

models_to_use <- c("q_constant", "b_constant", "b_rw", "b_ou")
###### setting #####

# data
alpha_increase_seq_1 <- seq(10, 750, by = 30)
alpha_decrease_seq_1 <- seq(750, 10, by = -30)
alpha_lamb =  c( rep(10,5), alpha_increase_seq_1 + rnorm(alpha_increase_seq_1,10,10), 
                 alpha_decrease_seq_1 + rnorm(alpha_decrease_seq_1,10,10),
                 rep(10,5))
beta_lamb = 0.5
T = 60
# reprot delay
D <- 15;

# Time period for checking
D_check <- 5

first_date <- as.Date("2024-01-01")

scoreRange <- c(first_date+days(9), first_date+days(19), first_date+days(29),
                first_date+days(39), first_date+days(49))

Fully reported case

Simulation

# q_constant_FR
params_q_constant_FR <- list(
  data = list(
    alpha_lamb = alpha_lamb,  
    beta_lamb  = beta_lamb,
    T       = T,
    date_start = first_date,
    D = D
  ),
  q_model = list(
    method        = "q_constant",
    method_params = list(b = 0.7, phi = 0.9)
  )
)

q_constant_FR <- simulateData(params_q_constant_FR)
plot(q_constant_FR$q, col = "red", pch = 19, type = "b", ylim = c(0,1))

Exploratory analysis

# exploritary analysis
page_num <- ceiling(nrow(q_constant_FR$case_reported_cumulated)/16)
exp_plot_q_constant <- fit_exp_plot(q_constant_FR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# print(exp_plot_q_constant)
# 
# exp_plot_q_constant$coefficients
exp_b_data_q_constant<- data.frame( date = as.Date(rownames(q_constant_FR$case_reported_cumulated)),
                          b = exp_plot_q_constant$coefficients$b)

exp_b_plot_q_constant <- ggplot(exp_b_data_q_constant, aes(x = date, y = b)) +
  geom_point(color = "black", size = 1.5) +       
  geom_smooth(method = "loess", se = TRUE,        
              color = "blue", fill = "grey", alpha = 0.5) +
  theme_minimal() +
  labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")

print(exp_b_plot_q_constant)
## `geom_smooth()` using formula = 'y ~ x'

Model fitting

# out_q_constant_FR <- nowcasting_moving_window(q_constant_FR$case_reported_cumulated, scoreRange =  scoreRange,
#                           case_true = q_constant_FR$case_true,
#                           start_date = first_date,
#                           D = D,
#                           methods =models_to_use,
#                           compiled_models = compiled_models,
#                           iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
#                           num_chains = 3, suppress_output = T,
#                           posterior_draws_path = file.path(posterior_draws_path, "q_constant")
#                           )
# 
# 
# save(out_q_constant_FR, file = file.path(data_save_path, "FR_q_constant.RData"))
load( file.path(data_save_path, "FR_q_constant.RData"))
# 
results_q_constant_FR <- nowcasts_table(out_q_constant_FR, D = D, report_unit = "day",
                          methods = models_to_use
                          )

plots_q_constant_FR <- nowcasts_plot(results_q_constant_FR, D = D, report_unit = "day",
                                  methods = models_to_use
                                  )
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
for (i in 1:length(results_q_constant_FR)) {
  print(calculate_metrics(results_q_constant_FR[[i]],
                          methods = models_to_use))
}
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 93.28 66.31 59.47 64.84         146.71           0.6 q_constant
## 2  6.99  9.26  4.35  6.52          64.61           1.0 b_constant
## 3  7.53  9.32  4.57  6.53          80.10           1.0       b_rw
## 4  8.20 10.20  4.95  7.16          75.51           1.0       b_ou
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 47.37 10.47 35.96 10.06         111.41           0.9 q_constant
## 2 18.34  2.95  8.61  2.15          87.40           1.0 b_constant
## 3 18.90  2.91  8.75  1.94          94.81           1.0       b_rw
## 4 19.40  2.93  8.38  1.90          93.16           1.0       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 22.75  2.74 16.15 2.54         113.54             1 q_constant
## 2 17.03  1.73  7.30 1.12         108.77             1 b_constant
## 3 19.65  1.81  7.73 1.20         118.37             1       b_rw
## 4 17.74  1.74  7.06 1.09         116.55             1       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 32.10  3.43 15.46 1.97         118.36             1 q_constant
## 2 25.47  2.78  7.96 1.17         115.94             1 b_constant
## 3 33.40  3.57  9.59 1.35         120.04             1       b_rw
## 4 25.10  2.77  7.71 1.13         117.81             1       b_ou
##    RMSE RMSPE  MAE MAPE Interval_Width Coverage_Rate     Method
## 1 15.03  3.09 9.59 1.66         110.35             1 q_constant
## 2 10.80  2.55 4.84 1.08         109.25             1 b_constant
## 3 11.57  2.79 4.64 1.08         109.91             1       b_rw
## 4 10.82  2.58 4.58 1.02         109.63             1       b_ou
for (i in 1:length(results_q_constant_FR)) {
  print(calculate_metrics(data.table::last(results_q_constant_FR[[i]], D_check),
                          methods = models_to_use))
}
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 131.14 65.07 104.67 62.36         257.63           0.4 q_constant
## 2   9.85 12.59   8.01  9.93         110.01           1.0 b_constant
## 3  10.63 12.81   8.53 10.25         140.60           1.0       b_rw
## 4  11.57 13.98   9.20 11.05         131.61           1.0       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 71.70 10.22 67.86 9.53         244.84           0.8 q_constant
## 2 36.27  4.55 27.82 3.66         191.01           1.0 b_constant
## 3 37.38  4.68 28.57 3.74         220.83           1.0       b_rw
## 4 38.50  4.79 28.12 3.67         214.43           1.0       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 39.79  2.77 33.77 2.38         267.82             1 q_constant
## 2 40.61  2.66 30.06 2.07         257.42             1 b_constant
## 3 47.27  3.06 32.71 2.22         311.02             1       b_rw
## 4 42.74  2.77 30.75 2.09         295.44             1       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 84.52  8.62 56.44 5.65         231.21             1 q_constant
## 2 71.37  7.29 45.73 4.57         224.63             1 b_constant
## 3 93.93  9.61 58.47 5.87         256.82             1       b_rw
## 4 70.49  7.20 45.12 4.51         239.02             1       b_ou
##    RMSE RMSPE   MAE MAPE Interval_Width Coverage_Rate     Method
## 1 38.27  8.77 28.29 6.29         147.62             1 q_constant
## 2 32.00  7.46 22.47 5.12         147.22             1 b_constant
## 3 35.05  8.23 22.94 5.34         153.61             1       b_rw
## 4 32.60  7.60 22.65 5.17         149.21             1       b_ou
plots_q_constant_FR
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

Parameter Checking

model q_constant

# try the third one, "2024-01-30"
T_test = T * 3/6

# q_constant
varnames_q_const <- out_q_constant_FR$q_constant[[3]]$summary()$variable

param_true = tibble(
    parameter = grep("^lambda\\[.+\\]$", varnames_q_const, value = TRUE),
    x = q_constant_FR$lambda[1:T_test]
)

mcmc_areas(out_q_constant_FR$q_constant[[3]]$draws("lambda"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^p\\[.+\\]$", varnames_q_const, value = TRUE),
    x = c(q_constant_FR$q[1], diff(q_constant_FR$q))
)

mcmc_areas(out_q_constant_FR$q_constant[[3]]$draws("p"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^q\\[.+\\]$", varnames_q_const, value = TRUE),
    x = q_constant_FR$q
)
mcmc_areas(out_q_constant_FR$q_constant[[3]]$draws("q"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^N\\[.+\\]$", varnames_q_const, value = TRUE),
    x = q_constant_FR$case_true[1:T_test, 1]
)
mcmc_areas(out_q_constant_FR$q_constant[[3]]$draws("N"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

model b_constant

# try the third one, "2024-01-30"
T_test = T * 3/6

# q_constant
varnames_b_const <- out_q_constant_FR$b_constant[[3]]$summary()$variable

mcmc_areas(out_q_constant_FR$b_constant[[3]]$draws("b"), prob_outer = 0.95)

mcmc_areas(out_q_constant_FR$b_constant[[3]]$draws("phi"), prob_outer = 0.95)

param_true = tibble(
    parameter = grep("^lambda\\[.+\\]$", varnames_b_const, value = TRUE),
    x =  q_constant_FR$lambda[1:T_test]
)

mcmc_areas(out_q_constant_FR$b_constant[[3]]$draws("lambda"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^q\\[.+\\]$", varnames_b_const, value = TRUE),
    x = q_constant_FR$q
)
mcmc_areas(out_q_constant_FR$b_constant[[3]]$draws("q"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^N\\[.+\\]$", varnames_b_const, value = TRUE),
    x = q_constant_FR$case_true[1:T_test, 1]
)
mcmc_areas(out_q_constant_FR$b_constant[[3]]$draws("N"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

Not fully reported

Simulation

# q_constant_NFR
params_q_constant_NFR <- list(
  data = list(
    alpha_lamb = alpha_lamb,  
    beta_lamb  = beta_lamb,
    T       = T,
    date_start = first_date,
    D = D
  ),
  q_model = list(
    method        = "q_constant",
    method_params = list(b = 0.1, phi = 0.9)
  )
)

q_constant_NFR <- simulateData(params_q_constant_NFR)
plot(q_constant_NFR$q, col = "red", pch = 19, type = "b", ylim = c(0,1))

Exploratory analysis

# exploritary analysis
page_num <- ceiling(nrow(q_constant_NFR$case_reported_cumulated)/16)
exp_plot_q_constant <- fit_exp_plot(q_constant_NFR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
# print(exp_plot_q_constant)
# 
# exp_plot_q_constant$coefficients
exp_b_data_q_constant<- data.frame( date = as.Date(rownames(q_constant_NFR$case_reported_cumulated)),
                          b = exp_plot_q_constant$coefficients$b)

exp_b_plot_q_constant <- ggplot(exp_b_data_q_constant, aes(x = date, y = b)) +
  geom_point(color = "black", size = 1.5) +       
  geom_smooth(method = "loess", se = TRUE,        
              color = "blue", fill = "grey", alpha = 0.5) +
  theme_minimal() +
  labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")

print(exp_b_plot_q_constant)
## `geom_smooth()` using formula = 'y ~ x'

Model fitting

# out_q_constant_NFR <- nowcasting_moving_window(q_constant_NFR$case_reported_cumulated, scoreRange = scoreRange,
#                           case_true = q_constant_NFR$case_true,
#                           start_date = first_date,
#                           D = D,
#                           methods =models_to_use,
#                           compiled_models = compiled_models,
#                           iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
#                           num_chains = 3, suppress_output = T,
#                           posterior_draws_path = file.path(posterior_draws_path, "q_constant")
#                           )
# 
# 
# save(out_q_constant_NFR, file = file.path(data_save_path, "NFR_q_constant.RData"))
load( file.path(data_save_path, "NFR_q_constant.RData"))

results_q_constant_NFR <- nowcasts_table(out_q_constant_NFR, D = D, report_unit = "day", 
                          methods = models_to_use
                          )

plots_q_constant_NFR <- nowcasts_plot(results_q_constant_NFR, D = D, report_unit = "day", 
                                  methods = models_to_use
                                  )

for (i in 1:length(results_q_constant_NFR)) {
  print(calculate_metrics(results_q_constant_NFR[[i]],
                          methods = models_to_use))
}
##    RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 24.31 24.80 18.22 22.68          94.21           1.0 q_constant
## 2 66.30 41.02 44.27 37.72          57.50           0.5 b_constant
## 3 72.26 43.68 48.16 40.15          61.00           0.5       b_rw
## 4 70.92 43.60 47.44 40.36          59.61           0.5       b_ou
##     RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 120.02 24.24 89.48 22.20          98.36          0.30 q_constant
## 2  76.18 16.46 54.25 15.28         118.32          0.75 b_constant
## 3  96.02 18.39 66.22 16.76         157.52          0.80       b_rw
## 4  82.14 17.00 58.28 15.83         146.62          0.85       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 187.91 21.64 144.55 20.47         108.40          0.20 q_constant
## 2  38.54 11.16  25.60  7.25         138.21          0.97 b_constant
## 3  44.87  9.44  29.19  7.00         243.13          1.00       b_rw
## 4  37.40 10.53  24.44  6.79         193.75          1.00       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 193.25 20.41 161.35 19.56         111.08          0.17 q_constant
## 2  37.82 10.28  26.93  6.32         137.36          1.00 b_constant
## 3  30.79  8.30  22.80  5.52         212.57          1.00       b_rw
## 4  34.59  9.96  23.56  5.86         178.71          1.00       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 185.44 20.75 157.32 19.94         102.59          0.16 q_constant
## 2  27.47  9.34  20.93  5.65         121.27          1.00 b_constant
## 3  34.51  8.85  24.17  5.85         155.83          0.98       b_rw
## 4  25.23  8.90  19.15  5.26         149.45          1.00       b_ou
for (i in 1:length(results_q_constant_NFR)) {
  print(calculate_metrics(data.table::last(results_q_constant_NFR[[i]],D_check),
                          methods = models_to_use))
}
##     RMSE RMSPE   MAE  MAPE Interval_Width Coverage_Rate     Method
## 1  33.89 18.54 31.16 18.31         162.42           1.0 q_constant
## 2  93.52 48.32 83.17 47.53          97.21           0.2 b_constant
## 3 101.96 52.51 90.85 51.84         104.01           0.2       b_rw
## 4 100.05 51.67 89.15 50.97         101.01           0.2       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 199.32 25.71 193.01 25.06         213.02           0.0 q_constant
## 2 128.61 16.60 117.95 15.28         258.24           0.8 b_constant
## 3 170.92 21.76 158.90 20.43         391.83           0.8       b_rw
## 4 140.98 18.11 130.80 16.90         343.64           0.8       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 350.03 24.54 344.44 24.28         251.01           0.0 q_constant
## 2  83.81  6.02  76.26  5.42         331.42           0.8 b_constant
## 3  93.44  6.67  86.99  6.17         705.27           1.0       b_rw
## 4  81.70  5.82  75.61  5.35         502.64           1.0       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 198.41 18.45 192.79 17.91         215.42           0.2 q_constant
## 2  79.50  7.73  71.94  6.84         276.82           1.0 b_constant
## 3  55.72  5.59  50.68  4.93         513.46           1.0       b_rw
## 4  71.23  6.72  59.57  5.61         393.62           1.0       b_ou
##     RMSE RMSPE    MAE  MAPE Interval_Width Coverage_Rate     Method
## 1 114.11 22.18 106.33 20.91         140.42           0.4 q_constant
## 2  43.82  9.58  37.29  7.92         170.21           1.0 b_constant
## 3  59.74 12.24  55.14 11.16         233.23           1.0       b_rw
## 4  41.26  8.94  35.22  7.42         201.01           1.0       b_ou
plots_q_constant_NFR
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

Parameter Checking

model q_constant

# try the third one, "2024-01-30"
T_test = T * 3/6

# q_constant
varnames_q_const <- out_q_constant_NFR$q_constant[[3]]$summary()$variable

param_true = tibble(
    parameter = grep("^lambda\\[.+\\]$", varnames_q_const, value = TRUE),
    x = q_constant_NFR$lambda[1:T_test]
)

mcmc_areas(out_q_constant_NFR$q_constant[[3]]$draws("lambda"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

param_true = tibble(
    parameter = grep("^p\\[.+\\]$", varnames_q_const, value = TRUE),
    x = c(q_constant_NFR$q[1], diff(q_constant_NFR$q))
)

mcmc_areas(out_q_constant_NFR$q_constant[[3]]$draws("p"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^q\\[.+\\]$", varnames_q_const, value = TRUE),
    x = q_constant_NFR$q
)
mcmc_areas(out_q_constant_NFR$q_constant[[3]]$draws("q"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^N\\[.+\\]$", varnames_q_const, value = TRUE),
    x = q_constant_NFR$case_true[1:T_test, 1]
)
mcmc_areas(out_q_constant_NFR$q_constant[[3]]$draws("N"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

model b_constant

# try the third one, "2024-01-30"
T_test = T * 3/6

# q_constant
varnames_b_const <- out_q_constant_NFR$b_constant[[3]]$summary()$variable

mcmc_areas(out_q_constant_NFR$b_constant[[3]]$draws("b"), prob_outer = 0.95)

mcmc_areas(out_q_constant_NFR$b_constant[[3]]$draws("phi"), prob_outer = 0.95)

param_true = tibble(
    parameter = grep("^lambda\\[.+\\]$", varnames_b_const, value = TRUE),
    x =  q_constant_NFR$lambda[1:T_test]
)

mcmc_areas(out_q_constant_NFR$b_constant[[3]]$draws("lambda"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^q\\[.+\\]$", varnames_b_const, value = TRUE),
    x = q_constant_NFR$q
)
mcmc_areas(out_q_constant_NFR$b_constant[[3]]$draws("q"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)

param_true = tibble(
    parameter = grep("^N\\[.+\\]$", varnames_b_const, value = TRUE),
    x = q_constant_FR$case_true[1:T_test, 1]
)
mcmc_areas(out_q_constant_NFR$b_constant[[3]]$draws("N"), prob_outer = 0.95) +
    geom_point(aes(x = x), param_true, color = "red", size = 1)